Michal Stawikowski Kohonen 1.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
from sklearn.metrics import accuracy_score
from sklearn.metrics import log_loss
import seaborn as sns
import warnings
from matplotlib import patches as patches
import copy
warnings.simplefilter("ignore")
%matplotlib inline
import math
plt.rcParams["figure.figsize"] = (12, 8)
plt.style.use("bmh")
font = {'family' : 'normal',
        'weight' : 'bold',
        'size'   : 22}

matplotlib.rc('font', **font)
from mpl_toolkits import mplot3d

Implementation

class Kohonen_Network(object):
    
    def __init__(self, M, N, data_dim, random = True):
        self.M = M #Width of map
        self.N = N #Height of map
        self.data_dim = data_dim #Input data dimension
        if(random):
            self.map = np.random.random((M, N, data_dim)) #Random map initialization
        else:
            self.map = np.zeros((M, N, data_dim))
        
    def lr_decay(self, i, epochs): #Learining decay
        return np.exp(-i / epochs)
    
    def gaussian(self, t, d): #Gausian function
        return np.exp(-((d*t*self.n_width)**2))
    
    def gaussian_s_der(self, t, d): #Gausin's function second derivative
        #return (2-4*(d*t*self.n_width)**2)*np.exp(-((d*t*self.n_width)**2))
        #return (2*((d*self.n_width)**2)-4*(((d*self.n_width)**2)*t)**2)*np.exp(-((d*t*self.n_width)**2))
        #return ((2-4*(d*self.n_width)**2)*np.exp(-((d*self.n_width)**2)))/2
        sigma = self.n_width*np.exp(-t / self.epochs)
        a = 1/((sigma**3)*math.sqrt(2*math.pi))
        b = (1-(d/sigma)**2)
        c = np.exp(-((d/sigma)**2)/2)
        return a*b*c
    
    def euclidean_metric(self,x, y):
        d = x - y
        return np.sqrt(np.sum(d * d, axis=-1))

    def closest_neuron(self, x):
        min_dist = np.inf
        min_indx = np.array([0,0])
        for i in range(self.M):
            for j in range(self.N):
                w = self.map[i,j]
                dist = self.euclidean_metric(w, x)
                if dist < min_dist:
                    min_dist = dist
                    min_indx = np.array([i,j])
        
        return (min_dist, min_indx)
    
    def train(self, data, epochs, n_fun = "gaussian", n_width = 1, lr = 1):
        self.n_width = n_width #Neighborhood size
        self.lr = lr #Learning rate
        self.epochs = epochs
        if(n_fun == "gaussian"):
            fun = self.gaussian
        else:
            fun = self.gaussian_s_der
            
        for t in range(epochs):
            shuffle = np.random.permutation(data.shape[0])
            data = data[shuffle,:]
            for x in data:
                min_dist, min_indx = self.closest_neuron(x)
                for i in range(self.M):
                    for j in range(self.N):
                        w = self.map[i,j]
                        dist = self.euclidean_metric(min_indx, [i,j])
                        self.map[i,j] = w + fun(t, dist)*self.lr_decay(t,epochs)*(x-w)
            
        
                    
    
    
    

Hexagon

data = pd.read_csv("C:/Users/jaiko/Desktop/MIOD/Kohonen/mio2/hexagon.csv",header=0)
#data.x = (data.x - data.x.mean())/data.x.std()
#data.y = (data.y - data.y.mean())/data.y.std()
x = np.array(data.iloc[:,0:2])
y = np.array((data.c))

Original data

matplotlib.rcParams.update(matplotlib.rcParamsDefault)
#sns.scatterplot(x = "x", y = "y", data = data, hue = "c")
plt.figure(figsize = [10,5]);
plt.scatter(x = "x", y = "y", data = data, c = "c");
plt.title("Data");
plt.show();

First try

net = Kohonen_Network(3,2,2)
net.train(x, epochs = 10, n_width = 1)
map_data = net.map
map_data = map_data.reshape(-1,2)

plt.figure(figsize = [10,5]);
plt.scatter(x = "x", y = "y", data = data, c = "c");
plt.scatter(map_data[:,0], map_data[:,1], c = "red", s = 100)
plt.title("Orginal data and 3x2 Kohonen map");
plt.show();

As we can see it detected clusters pretty well.

Gaussian function

neigh_size = np.arange(0.1,1.1,0.3)
font = {'family': 'serif',
            'color':  'black',
            'weight': 'normal',
            'size': 10,
            }
for size in neigh_size:
   
    net = Kohonen_Network(10,10,2)
    net.train(x, epochs = 10, n_width = size)
    map_data = copy.copy(net.map)
    map_data = map_data.reshape(-1,2)

    plt.figure(figsize = [10,5]);
    plt.scatter(x = "x", y = "y", data = data, c = "c",s = 80, alpha = 0.5);
    plt.scatter(map_data[:,0], map_data[:,1], c = "red");
    i=0
    for a in range(1, net.M + 1):
        for b in range(1, net.N + 1):
            plt.text(net.map[a-1,b-1,:][0], net.map[a-1,b-1,:][1], str(i), fontdict=font)
            i += 1
    plt.title("Neighborhood size: " + np.array2string(size));
    
    plt.show();
pass;   

As can be seen from the above charts, clusters have more or less been identified, their number corresponds to the number of classes, and their location is similiar to these from the original data. With the growing size of the neighborhood, the uneven density of points in individual clusters increased. Network learning time is quite long compared to MLP. Neurons close at MxN map have simmiliar weights as we can see from indexing.

Weight Heatmap

x = np.array(data.iloc[:,0:2])
y = np.array((data.c))
net = Kohonen_Network(10,10,2)
net.train(x, epochs = 3, n_width = 0.4)

map_data = copy.copy(net.map)
map_data = map_data.reshape(-1,2)
map_data[:,0] = (map_data[:,0] - np.min(map_data[:,0]))/(np.max(map_data[:,0])-np.min(map_data[:,0]))
map_data[:,1] = (map_data[:,1] - np.min(map_data[:,1]))/(np.max(map_data[:,1])-np.min(map_data[:,1]))

col = map_data.reshape(net.M,net.N,2)
fig = plt.figure()
zeros = np.zeros((net.M,net.N,3))
zeros[:,:,2] = 0.5
zeros[:,:,:2] = col


ax = fig.add_subplot(111, aspect='equal')
ax.set_xlim((0, net.M+1))
ax.set_ylim((0, net.N+1))
ax.set_title('Kohonen Map')

# plot
for x in range(1, net.M + 1):
    for y in range(1, net.N + 1):
        ax.add_patch(patches.Rectangle((x-0.5, y-0.5), 1, 1,
                     facecolor=zeros[x-1,y-1,:],
                     edgecolor='none'))
plt.show()

We can clarly see the division into all clusters here and that close neurons have simlilliar weights

Gaussian second derivative

data = pd.read_csv("C:/Users/jaiko/Desktop/MIOD/Kohonen/mio2/hexagon.csv",header=0)
data.x = (data.x - data.x.mean())/data.x.std()
data.y = (data.y - data.y.mean())/data.y.std()
x = np.array(data.iloc[:,0:2])
y = np.array((data.c))

net = Kohonen_Network(10,10,2)
net.train(x, epochs = 8, n_width = 5, n_fun = "second")

map_data = copy.copy(net.map)
map_data = map_data.reshape(-1,2)
map_data[:,0] = (map_data[:,0] - np.min(map_data[:,0]))/(np.max(map_data[:,0])-np.min(map_data[:,0]))
map_data[:,1] = (map_data[:,1] - np.min(map_data[:,1]))/(np.max(map_data[:,1])-np.min(map_data[:,1]))

col = map_data.reshape(net.M,net.N,2)
fig = plt.figure()
zeros = np.zeros((net.M,net.N,3))
zeros[:,:,2] = 0.5
zeros[:,:,:2] = col


ax = fig.add_subplot(111, aspect='equal')
ax.set_xlim((0, net.M+1))
ax.set_ylim((0, net.N+1))
ax.set_title('Kohonen Map')

# plot
for x in range(1, net.M + 1):
    for y in range(1, net.N + 1):
        ax.add_patch(patches.Rectangle((x-0.5, y-0.5), 1, 1,
                     facecolor=zeros[x-1,y-1,:],
                     edgecolor='none'))
plt.show()   

Individual clusters are separated from each other by negative values, which is produced by the second derivative of the Gaussian function.

Cube

data = pd.read_csv("C:/Users/jaiko/Desktop/MIOD/Kohonen/mio2/cube.csv",header=0)
#data.x = (data.x - data.x.mean())/data.x.std()
#data.y = (data.y - data.y.mean())/data.y.std()
#data.z = (data.z - data.z.mean())/data.z.std()
x = np.array(data.iloc[:,0:3])
y = np.array((data.c))

Original data

plt.figure(figsize=(20,10))    
ax = plt.axes(projection='3d')
ax.scatter(data.x, data.y, data.z, c=data.c)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

Text(0.5, 0, 'z')

First try

net = Kohonen_Network(4,2,3)
net.train(x, epochs = 20, n_width = 0.1)
map_data = net.map
map_data = map_data.reshape(-1,3)

plt.figure(figsize = [10,5]);
ax = plt.axes(projection='3d')
ax.scatter(data.x, data.y, data.z, c=data.c, alpha = 0.1)
ax.scatter(map_data[:,0], map_data[:,1], map_data[:,2], c = "red", s = 100)
plt.title("Orginal data and 4x2 Kohonen map");
plt.show();

Gaussian function

font = {'family': 'serif',
            'color':  'black',
            'weight': 'normal',
            'size': 20,
            }
neigh_size = np.arange(0.1,1.1,0.3)

for size in neigh_size:
    i=0
    net = Kohonen_Network(7,7,3)
    net.train(x, epochs = 10, n_width = size)
    map_data = net.map
    map_data = map_data.reshape(-1,3)

    plt.figure(figsize=(30,20)) 
    ax = plt.axes(projection='3d')
    ax.scatter(data.x, data.y, data.z, c=data.c, alpha = 0.3)
    ax.scatter(map_data[:,0], map_data[:,1], map_data[:,2], c = "red", s = 100)
    for a in range(1, net.M + 1):
        for b in range(1, net.N + 1):
            ax.text(net.map[a-1,b-1,:][0], net.map[a-1,b-1,:][1], net.map[a-1,b-1,:][2], str(i), fontdict=font)
            i += 1
    ax.set_xlabel('x', fontsize = 40)
    ax.set_ylabel('y', fontsize = 40)
    ax.set_zlabel('z', fontsize = 40)
    plt.title("Neighborhood size: " + np.array2string(size), fontsize = 40);
    plt.show();
pass;  

In three dimensions it is harder to see, but the network has divided the points in a similar way as they appear in the original set. The higher the size of the neighborhood, the more widespread points are.

Weight Heatmap

data = pd.read_csv("C:/Users/jaiko/Desktop/MIOD/Kohonen/mio2/cube.csv",header=0)
x = np.array(data.iloc[:,0:3])
y = np.array((data.c))
net = Kohonen_Network(20,20,3)
net.train(x, epochs = 4, n_width = 0.3)
map_data = net.map
map_data = map_data.reshape(-1,3)
map_data[:,0] = (map_data[:,0] - np.min(map_data[:,0]))/(np.max(map_data[:,0])-np.min(map_data[:,0]))
map_data[:,1] = (map_data[:,1] - np.min(map_data[:,1]))/(np.max(map_data[:,1])-np.min(map_data[:,1]))
map_data[:,2] = (map_data[:,2] - np.min(map_data[:,2]))/(np.max(map_data[:,2])-np.min(map_data[:,2]))
fig = plt.figure()
col = map_data.reshape(net.M,net.N,3)
ax = fig.add_subplot(111, aspect='equal')
ax.set_xlim((0, net.M+1))
ax.set_ylim((0, net.N+1))
ax.set_title('Kohonen Map')

# plot
for x in range(1, net.M + 1):
    for y in range(1, net.N + 1):
        ax.add_patch(patches.Rectangle((x-0.5, y-0.5), 1, 1,
                     facecolor=col[x-1,y-1,:],
                     edgecolor='none'))
plt.show()

We can see that there is happening a lot. Nuerons that are close together have similar weights, you can see a division into several clusters.

Gaussian second derivative

data = pd.read_csv("C:/Users/jaiko/Desktop/MIOD/Kohonen/mio2/cube.csv",header=0)
data.x = (data.x - data.x.mean())/data.x.std()
data.y = (data.y - data.y.mean())/data.y.std()
data.z = (data.z - data.z.mean())/data.z.std()
x = np.array(data.iloc[:,0:3])
y = np.array((data.c))
net = Kohonen_Network(20,20,3)
net.train(x, epochs = 4, n_width = 7, n_fun = "second")
map_data = net.map
map_data = map_data.reshape(-1,3)
map_data[:,0] = (map_data[:,0] - np.min(map_data[:,0]))/(np.max(map_data[:,0])-np.min(map_data[:,0]))
map_data[:,1] = (map_data[:,1] - np.min(map_data[:,1]))/(np.max(map_data[:,1])-np.min(map_data[:,1]))
map_data[:,2] = (map_data[:,2] - np.min(map_data[:,2]))/(np.max(map_data[:,2])-np.min(map_data[:,2]))
fig = plt.figure()
col = map_data.reshape(net.M,net.N,3)
ax = fig.add_subplot(111, aspect='equal')
ax.set_xlim((0, net.M+1))
ax.set_ylim((0, net.N+1))
ax.set_title('Kohonen Map')

# plot
for x in range(1, net.M + 1):
    for y in range(1, net.N + 1):
        ax.add_patch(patches.Rectangle((x-0.5, y-0.5), 1, 1,
                     facecolor=col[x-1,y-1,:],
                     edgecolor='none'))
plt.show()

Individual clusters are separated from each other by negative values, which is produced by the second derivative of the Gaussian function.

Źródła:

• http://mini.pw.edu.pl/~karwowskij/mioad/wyklad2.pdf - lectures
• https://github.com/tcosmo/tcosmo.github.io/blob/master/assets/soms/code/SOM_viz.ipynb - idea for visualisation
• https://books.google.pl/books?id=-fQGKI8ke_cC&pg=PA76&lpg=PA76&dq=som+second+derivative+gaussian&source=bl&ots=mIXex-tGMG&sig=ACfU3U2tHFMHlOWq28loSzx0-kFCEXDktA&hl=en&sa=X&ved=2ahUKEwjDm5WsvonpAhVKkMMKHSD8A20Q6AEwCHoECAwQAQ#v=snippet&q=sigma&f=false - second derivative idea

,,Potwierdzam samodzielność powyższej pracy oraz niekorzystanie przeze mnie z niedozwolonych źródeł"